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CHAPTER 1: INTRODUCTION 



The puq)ose of this thesis is to study the stability characteristics of a pair of coupled van 
der Pol equations. In order to fulfill this purpose, several different concepts must be 
discussed prior to actually looking at the behavior of the coupled equations. 

The first issue to be dealt with is the van der Pol equation itself; this topic will be covered 
in this chapter with details of the solution to the van der Pol equation via perturbation 
theory given in Chapter 2. 

Chapter 3 will explore the derivation, existence and computation of Pade approximants 
which will be used to model behaviors of the van der Pol equation and the coupled 
equations. If the reader is already comfortable with the characteristics of Pade 
approximants. Chapter 3 can be skipped. 

In Chapter 4, the model for the coupled van der Pol oscillators will be examined. The 
scope of the problem of characterizing the stability behaviors of the coupled equations will 
be examined. To that end, the Zero Mean Damping (ZMD) surface will be introduced and 
the variational equations to the coupled oscillators will be derived. Some implications of 
Floquet theory on the variational equations will also be explored. 

Chapter 5 will detail the solution process of the variational equation derived in Chapter 4. 
The stability transition curves along the ZMD surface will be introduced and the method 
for determining them will be presented. 

Chapter 6 will summarize the results of Chapter 5 and discuss some ongoing work and 
computational difficulties encountered in determining the stability transition curves. 
Further studies will also be discussed. 



TFE VAN DER POL EQUATION 
The van der Pol equation: 



u-e(l-u^)u + u = 0 



(1-1) 
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was introduced by Bismarck van der Pol in a paper published in 1927 [29,30] to model the 
non-linear resistance of certain circuits including vacuum tubes. Since that time Equation 
(1-1) has become a popular mathematical model for limit cycle oscillators 

A limit cycle oscillator, is an oscillator which exhibits self-sustained periodic behavior. It 
is non-conservative, due to a damping term, and non-linear [2,7]. In the case of the van 
der Pol equation, the non-linearity is also due to the damping term. Limit cycle behavior is 
present in many mechanical, electrical, and biological systems [3 ,4, 5, 6]. Stoker[21] and 
LaSalle[15] showed the existence of a unique limit cycle for Equation (1-1) for all e>0. 
The limit cycle in the phase plane (u, u’} is shown in Figure (1-1) for the case of 8=0.9. 

The shape of the limit cycle for the van der Pol equation varies markedly with the value of 
the non-linear parameter, 8. For 8«1 the limit cycle is well approximated by u=2 Cos(t), 
while for large 8 the limit cycle remains periodic but varies significantly from sinusoidal 
[10,28,29]. For large 8, the van der Pol equation can be used to model relaxation 
oscillators such as a heart beating, a clock ticking or a fish swimming [6,1 1,22,28]. 




Figure (1-1) Phase-plane for VDP 
Oscillator (8=0.9). 

Applications of relaxation oscillators is the motivation for this thesis. If the behavior of a 
chain of linked van der Pol equations can be characterized, it may be possible to build an 
electro-mechanical system which uses coupled oscillators to produce control signals for 
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locomotion. As an example, consider a pair van der Pol oscillators in the relaxation limit. 
Since u(t) is characterized by a rapid transition between two relatively flat states, using the 
signal to drive a pair of servomechanisms may be a very inexpensive way to drive a bipedal 
robot. By coupling the oscillators, it may be possible to produce a smooth, responsive 
gait with little or no computer processing. 

Another possible application for coupled van der Pol oscillators could be modeling the 
synapses in a fish’s spine while swimming [6]. Here the chain of oscillators may be 30 or 
even 50 units long and modeling may be extremely complicated. 



In the first step to understanding the behavior of such models, it is desirable to find the 
stability characteristics of the coupled equations. In this thesis, the scope will be limited to 
only two oscillators which could be modeled as; 



X - s(l - x^ )x + X = sA(y - x) + eB(y - x) 
y - s{l - y^ )y + y = sA(x - y) + sB(x - y) 



(1-2) 



where sA represents the displacement coupling and sB represents the velocity coupling 
between the two oscillators. This model was proposed by Rand and Holmes [18] and 
further studied by Storti and Reinhall[ 19,20,22,27], In future, examining the stability of a 
chain of more than two oscillators could be done by extending the techniques of Storti and 
Reinhall and this thesis. 



CHARACTERIZING THE VAN DER POL EQUATION 

The first step in understanding the stability of Equations (1-2) is to develop the solution to 
the van der Pol equation. Equation (1-1) has been studied extensively [1, 2, 7, 8] and a 
simple approach to finding solutions using Lindstedt’s method [14] is presented in Chapter 
2 . 



^ A ■ 




CHAPTER 2: SOLUTION TO THE VAN DER POL EQUATION 



Perturbation theory is one of the most powerful techniques used in solving non-linear 
differential equations in vibration theory [14,17], In concept, perturbation theory is very 
simple and can be explored by examining the van der Pol equation in time t: 

u-8(l-u^)u + u = 0, u(0) = 0 (2-1). 

Using Lindstedt’s method, rescale time and let x = o t. Changing variables. Equation (2- 
1) becomes: 



u"-s(l-u')ou' + o)'u = 0, u'(0) = 0 (2-2) 

where the primes represent differentiation with respect to the rescaled time t. Now apply 
perturbation theory, by considering Equation (2-2) as a small perturbation of the linear 
system; 



u"+u = 0, u'(0) = 0 (2-3) 

whose behavior is well understood. It is reasonable to assume that u(t) and co may be 
represented by power series in e; 

U(t) = Uo(t) + 8U,(t) + 8^U2(t)+--- 

ft) = 1+8 (D, +8^(02 +••• (2-4) 

Substituting Equations (2-4) into Equation (2-2) yields: 

Uo +8 u,"+8^u" + •••+(! + 8 ft), + 8^(0 2 +*'*)^(Uo +8 U, +8^U2 +•••) = 

8 (l -(Uq +8 U, +8^U2 H )^}(1 +8 ft), + 8 ^ft )2 + •••)(uj +8 Uj +8^0^ H ) 

(2-5) 



Collecting like powers of 8, the first few equations become: 
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8«; u" + u„=0 

e': u"+u, =u;(l-Uo^)-2w,Uo 

u" + U2 =co,u;(l-Uo^) + u;(l-Uo^)-2u,UoU; -2(0, u, -®,^Uo 



(2-6). 



These equations are then solved, in order, where the right hand side is completely fixed by 
the solutions to the previous equations and determining the (Oj’s and any undetermined 
coefficients in the earlier u i(t)’s to suppress secular terms. 

Secular terms are terms that if not canceled in some manner, would give solutions to the 
differential equation that grow with time. As an example, consider the linear Equation (2- 
4) Avith an added driving force; 

u" + u = Cos(t) (2-7). 

Since the driving force matches the natural frequency, u(t) will grow over time and 
become unbounded; the solution to Equation (2-7) is u(t)=l/2 1 sin(t), which clearly grows 
over time. In applying perturbation theory, the goal is to select the undetermined 
coefficients in Equation (2-4) to eliminate terms which match the frequency of Equation 
(2-3). In particular, for the van der Pol equation, the coefficients of any Sin(t) or Cos(t) 
term must be suppressed when they appear on the right hand sides of Equations (2-6). 

As a demonstration, consider the first few orders of 8 in Equations (2-6). Starting with 
zeroeth order, the solution must be of the form: 

Uo(t) = AoCos(t) + BoSin(t) (2-8). 

Now applying the initial condition in Equation (7), Bo=0 and Ao will be determined while 
solving the first order equation. Substituting xo(t) into the first order equation and 
employing trigonometric identities; 

uj'+u, = (AoCos(t))'(l“(AoCos(t))^)-2cc>,(AoCos(t)) 

f 1 ^ 1 ( 2 ’ 9 ). 

= |^-Ao+ — AoJ Sin(t) + — Aq Sin(3t)-2co,AoCos(t) 

In order to suppress the secular terms, the coefficients of Sin(t) and Cos(t) must be set 
equal to zero; there are three possible values for Ao to do this, Ao=-2, 0 or 2 and ( 0 i must 
be zero. Selecting Ao=0 will result in the trivial solution for the zeroeth order equation. 
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but either 2 or -2 would be valid and give solutions to Equation (2-1) with opposite signs. 
For convenience, take the solution Ao=2. The right hand side of Equation (2-9)is now 2 
Sin(3t) and the solution is; 

u,(t) = B,Sin(t)-^Sin(3t) (2-10) 

and imposing the initial condition of Equation (2-1) again: 

u,(t) = ;^Sin(t)-^Sin(3t). (2-11) 

It serves little point to continue this procedure here, but it is worth noting that both the 
first order term of the solution, ui(t), and the first order term in the frequency expansion, 
toi, were determined from the first order equation of Equations (2-6). This pattern will 
hold while solving Equations (2-6) at each order of s. Following this procedure, it is 
possible to find solutions for all Uj(t)and coj up to a high value of i; such a solution is 
considered an i“ order solution to the differential equation [14]. 

This can be a tedious solution procedure, but the main advantage is that any degree of 
accuracy can be achieved. In particular, if the solution is needed to some high order n, 
this method will allow determining the solution of Equation (2-1) through order n with 
sufficient time and effort. Fortunately, this method (Lindstedt’s method) can be realized 
using a programmable, symbolic algebra package. The results outlined here were 
determined using a program written in Mathematica™ [9, 32]; the code to generate these 
results is found in Appendix A. This code takes advantage of some additional features of 
the van der Pol equation that will be discussed in the next section. 

SYMMETRY IN THE VAN DER POL EQUATION 

Examining Equation (2-6) with symmetry of the solution in mind, it can be seen that the 
right hand sides of the equations alternate in symmetry. In particular, the right hand side 
of the zeroeth order equation is 0 which is even, therefore uo(t) must be even, which it is. 
The right hand side of the first order equation, considering the fact that Uo(t) is even, must 
be odd resulting in an odd solution for ui(t). This pattern continues and allows the rapid 
solving of each order of the power series solution. 

One additional property of the equation may be used to accelerate the solution of 
Equations (2-6). Since the left hand side of every order in s is the same, namely; 



LHS = u"+Uj 



(2-12) 
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then the homogeneous solution is always of the form: 

u,,(t) = A, Cos(t) + BiSin(t) (2-13) 

and due to the argument of the above paragraph, Aj=0 if i is odd or Bi=0 if i is even. 

Even more fortunately, the particular solution for each order in e is also predictable. The 
right hand side of Equation(2-6), after suppressing secularity, will always take the form 
of 



(2i+l),A=2 

RHS = ao + X Cos(k t) Vie Even 

k=3 

(2i+l),A=2 

RHS= 2bkSin(kt) VieOdd 

k=3 

This results in particular solutions of the form: 

(2i+l),A=2 ^ 

H 2 “ ,Cos(kt) VieEven 

k=3 IC “ I 
(2i+l),A=2 

’^i.p(^)= S VieOdd 

k=3 IC — I 



(2-14). 



(2-15). 



The overall advantage of being able to write down these solutions, is that Lindstedt’s 
method may now be used without actually solving any differential equations. This 
significantly speeds up the solution generation for the van der Pol equation. 

After removing the need for using a differential equation solver, the next most time- 
intensive process left in solving each order of the van der Pol equation is expanding the 
right hand side of each equation and using the trigonometric identities to reduce them to 
the form of Equation (2-14). One very simple way to reduce the computational effort is to 
use Euler’s identity; 



e^ ‘ = Cos(t) + j Sin(t), or 
e^‘ +e“^‘ 

Cos(t) = r and Sin(t) = 



e'* +e-j‘ 



2 



2j 



(2-16) 
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where j = . Further, in using Mathematical^, the expression Exp(j t) is very 

complicated, so it is advantageous to define another symbol, say Z = Exp(j t), and then 
replace each sine and cosine term on the right hand side of Equations (2-6) with: 



Cos(kt) 



+Z‘‘ 

2 



and 



Sin(k t) 




(2-17) 



Using these expressions, the right hand side equations can be expanded much more rapidly 
and then regrouped back into terms of sines and cosines to implement Equation (2-15) and 
suppress secularity. 



The first few terms of the solution of Equation (2-1) are presented here: 

1 \ 2 f Cos (t) 



u (t) = 2 Cos (t) + 
/ 3 Sin (t) 

n— 



3 5 ^ 

+ — Cos (3 1) Cos (5 1)| + 

16 Q6 V 



,, 7Sin(t) 21 
( tt^+ — Sin( 



256 



256 



35 7 

in (3 1) Sin (5 1) + Sin (7 1)) + O (e"*) 

576 576 > 



(2-18) 



with the frequency expansion given by: 

17€“ 35e^ 678899e* 28160413e'<^ , ,, 

w = 1 + + + + 0(e'^ 

16 3072 884736 5096079360 2293235712000 ^ ’ (2-19). 

Using the short-cuts listed above, the code found in Appendix A was developed. The first 
ten terms in the solution of the van der Pol equation are found in Appendix C along with 
numerical solutions up through order 25. The solution to the van der Pol equation was 
determined analytically (no numerical rounding) up through order 75 using this code. 



CHAPTER 3: BRIEF DESCRIPTION OF PADE APPROXIMANTS 



WHAT IS A PADE APPROXIMANT? 

A Pade approximant is a rational function which approximates a power series [12,13], It 
is represented by [L/P] where the numerator is of degree L and the denominator is of 
degree P. The Pade approximant [L/P] can be constructed from any power series of 
degree N; 



f(z) = ^CjZ' = Cq + C,Z + C 2 Z^+-*-+Cn z^ + o(z^^’) 



(3-1) 



where L + P < N, and takes the form; 



[L / P] = ^o+a.z + a,z^+-;-H^z ^ j 

b„ + b, z + b 2 z +• • -i-bp z 



(3-2) 



and the Pade approximant’s Maclaurin series expansion must match Equation (3-1) 
exactly up through order L + P. In order to make [L/P] unique (since multiplying the 
numerator and denominator by any value would create another approximant with the same 
Maclaurin series), it is common to arbitrarily set bo=l (or divide both the numerator and 
denominator by bo), thus making Equation (3-2): 



[L/P] = 



a<) +a,z + a2Z^+-"+aLZ^ 
1 + b,z + b2Z^+-”+bpZ^ 



+ 0(z™) 



(3-3) 



Note: The power series of Equation (3-1) is assumed to be complex- 

valued and could represent a Maclaurin series expansion of some function 
f(z); however, the power series need not be complex-valued in order to 
work with Pade approximants and Equation (3-1) could just as easily 
represent the Taylor’s series expansion of a real valued function, or an 
experimentally obtained polynomial curve fit, or a power series solution 
obtained with perturbation theory. 



WHAT ARE FADE APPROXIMANTS FOR? 

There are three main reasons for constructing Fade approximants. 



10 



The first reason to construct [L/P] is to produce a function which has a larger radius of 
convergence than the original power series. In some cases [L/P] may have a significantly 
larger radius of convergence than the original power series, or possibly even converge for 
all values of z. Even if the approximant does not have a larger radius of convergence, the 
approximant may yield information about how the function behaves for large values of z. 
This will be more readily seen in an example. 

Example 3-1 

Consider a power series: 

f(x)= 1-x + x^ - x^ + ••• 

which is the Taylor’s series for 1/(1 + x) has a radius of convergence R<1. 

But the an approximant for f(x) is given by 



which is the original function and is finite for all x except x=-land 
converges for R<oo. One can see the behavior of f(x) and [1/1] easily in a 
graph (f(x) is the dotted line, [1/1] is solid) where f(x) includes the first 8 
terms in the series: 

f{x) & [1/1] 




Taylor’s Series 



3 



2 




0 



X 



1 



2 



3 



Figure (3-1) Taylor’s series and Fade approximant for l/(l+x). 
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Here, the approximant obviously does a better job of estimating the 
function for x>l, since the approximate is exactly the same as the original 
function. 



The second reason for using Fade approximants is to speed up the evaluation of a given 
power series when working inside the radius of convergence. For the preceding example, 
it is less computationally expensive to find the value of 1+ x and its reciprocal than it is to 
evaluate the Taylor’s series. In a more general sense, comparing the evaluation of two 
polynomials of degree N/2 and taking their ratio is a less costly evaluation than calculating 
the value of an N degree polynomial, due to the expense of calculating higher powers in 
the N degree polynomial. 

The third reason for employing Fade approximants, which is closely related to the first, is 
to more accurately model the asymptotic behavior of a function. One of the worst 
features about a power series solution is that the highest order term will always dominate 
when X is large, so the series behaves as x’^. This may be unfortunate, especially if there is 
some indication of how the actual solution should behave. For example, if the function 
being modeled is known to decay as 1/x or is known to approach a constant value, then a 
Fade approximant should be more successful in predicting the asymptotic behavior by 
careful selection of L and F when constructing the [L/F] approximant. In the case of a 
function which should decay as 1/x, selecting F=L+1 will probably result in an 
approximant that also behaves as 1/x in the limit as x becomes large since the numerator 
will behave as x'' and the denominator as x^^V In the case where the function should 
approach a constant, picking L=F will most likely produce an approximant with the 
desired behavior. (Example 3 below will explore some of this more clearly.) 

CALCULATING FADE AFFROXIMANTS 



Frocedure 



Calculating a Fade approximant is a straight-forward process [13]. Starting by setting 
Equation (3-1) equal to Equation (3-3); 



Co + C, Z + C^ z^ +• • -t-CN z^ 



ao +a,z + a2Z^-i---+aLz‘' 
l + bjZ + b2Z^-l“' •+bp z^ 



Next, multiply both sides by the denominator of the right side; 



(3-4) 



(co + c, z -I- C 2 z^ 4- • -f-Cj^ z’^)(l -I- b, z + b 2 z^ 4- • -4-bp z**) = 



ao 4-a,z4-a2Z^4-'--4-aLZ^ 



(3-5) 



12 



Multiplying Equation (3-5) out and collecting like powers of z, yields: 



^L-P+1 ^L-P+2 




'bp “ 






^L-P+2 ^L-P+3 ^L+1 




bp_i 




^L+2 


: 




. 


= — 


. 


• i \ i 




• 




• 


^L-P+1 ^L+P-l_ 








_^L+P_ 



where cj = 0 ifj<0 for consistency, and 



a, =c, +b,c„ 

= c^ +b,c, H-bjCo 

min(L,P) 

ZbiCL_i 



(3-6) 



(3-7). 



The linear system in Equation (3-6) may be solved for the b;’s by any suitable method, and 
then the bi’s used in Equations (3-7) to find the ai’s. Following this procedure, any order 
Fade approximant can be constructed for L + P < N. 

Equations (3-6) and (3-7) can be easily coded into any programming language or written 
as functions in many mathematics packages. 



EXAMPLES 



Example 3-2 

Find [4/4] for the Taylor’s series approximation to Sin(x). Starting with 
the Taylor’s series: 



f(x) = x-- + --- + 0(x ) 



and [4/4] is 



[4/4] 



a,) +a,x + a2X^ -i-a^x'* 

l + b,x + b2x^ -i-bjX^ +b4x'* 
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Using Equation (3-6): 



1 


0 


-1/3! 


0 


X' 




■ 1/5! ■ 




’^ 4 ' 




‘11/5880' 


0 


-1/3! 


0 


1/5! 


bj 




0 




t>3 




0 


1/3! 


0 


1/5! 


0 


b. 


— “ 


-1/7! 


=> 


b. 


— 


3/49 


0 


1/5! 


0 


-l/7!_ 


b, 




0 




.b|. 




0 



a« =0 

a, =1 + 0 = 1 

2i2 = 0 + 0 + 0 = 0 

aj = -1/6 + 0 + 3/49 + 0 = -31/294 
a4 = 0 + 0 + 0 + 0 + 0 = 0 

so the final solution is: 



X- 



[4/4] = 



31 

294 



1 + — + 
49 



11 

5880^ 



O(x’) 



Note: Fade approximants do not perform much better than regular power 
series in predicting periodic Sanctions; this example was included for 
numerical interest as will be discussed later. 

Example 3-3 

Another example of finding a Fade approximant demonstrates how 
asymptotic behavior for Fade approximants may be significantly better than 
a normal power series. Consider the fianction: 



f(x) = with Fade approximant [3 /4] = 1 — . 

+ 1 + x^+V 

8 

Selecting a [3/4] approximant increases the likelihood of correctly 
modeling the 1/x behavior of f(x) for large x. Attempting to calculate the 
[3/4] approximant yields a [2/4] approximant which shows that not all 
degree approximants are necessarily unique. Below is a graph of f(x) (solid 



line), [3/4] (dotted line) and the eighth order Taylor’s series for f(x) (dash- 
dot). In this case, the [3/4] approximant does a significantly better job of 
modeling the original equation and tends to zero (as 1/x^) for large x. 




Figure (3-2) Taylor’s series and a [3/4] approximant to f(x) = ■ , ' • 

Vl + x 



STABILITY OF FADE APPROXIMANTS 
Condition Number Considerations 

Following the procedure set out here for finding a Fade approximant includes solving the 
linear system represented by Equation (3-7) and Equations (3-7). While Equations (3-7) 
present no difficulties computationally, the solution to Equation (3-6) may be subject to 
significant loss of numerical accuracy. This is of particular concern when combined with 
the fact that the accuracy of Fade approximants is highly sensitive to the differences 
between coefficients in the approximant. (The need for very high numerical accuracy is 
discussed in detail in Reference [13], section 2.1.) 

To illustrate how significant the numerical results may be affected when solving Equation 
(3-6), one can examine the condition number (or a reasonable approximation for the 
condition number) of the matrix C: 
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^L-P+l ^L-P+2 
_ ^L ^L+P+1 

To gain some measure of the magnitude of the condition number of this matrix, consider 
the calculations carried out in Example 3-2. Below is a table which shows how the 
condition number of C varies with the degree of the Fade approximant used to estimate 
the function f(x)=Sin(x) and the number of lost digits in accuracy is assumed to be 
Logio(condition number): 



'L+i 



'L+P-1 



(3-8), 



L=P Cond # Lost Digits 



4 


4.72 


lO" 


4 


5 


1.66 


10^ 


5 


6 


4.27 


10’ 


8 


7 


3.40 


10^ 


10 


8 


1.74 


lO’’ 


12 


9 


2.46 


lO*'* 


14 


10 


2.08 


lO” 
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Table (3-1) Condition number and loss of 
accuracy for Example (3-2). 



The importance of finding an accurate method to solve Equation (3-6) is apparent with 
respect to these results, and it would be wise to check the condition number of C when 
employing Fade approximants. 

There is one obvious way to avoid the numerical problems discussed here. By employing 
one of the powerful math packages available which can do exact arithmetic such as 
Mathematica™ or Maple™, Equation (3-6) can be solved with no loss of accuracy when 
the coefficients cj are known as precise fractions. While it may seem optimistic to want to 
know the cj’s exactly, in many applications where Fade approximants are used, the power 
series being replaced is generated exactly from a technique such as perturbation theory. In 
particular, in non-linear vibration analysis, the solution to a differential equation may be a 
power series in the non-linear parameter, 8, known to 50 or even 100 terms which are 
calculated with fractions and known precisely. This will be discussed more fully when 
examining applications of Fade approximants to vibration analysis. 



What is a defect? 
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In simple terms, a defect is a failure of a Fade approximant to accurately approximate a 
power series in the neighborhood of a root of the denominator of the approximant. Or in 
mathematical terms, a defect occurs around x==xd when l+b,XD +b 2 Xp+---+bpXD =0 . 
Borrowing on the language of complex analysis, the defect occurs where [L/P] has a pole 
at xd or zd. 

How to handle a defect 

The presence of the pole in [L/P] may or may not be a result of numerical inaccuracy in 
calculating the approximant, and on a first look it seems reasonable that there would be 
cases where we would expect the poles to be present. The most obvious reason for the 
[L/P] approximant to have a pole would be if the original function, f(z), also has a pole. 
In this case, the behavior of the approximant will depend on the nature of the pole in f(z); 
if the pole is an essential pole (there is no finite n such that (z-zd)" f(z) is analytic at zd) 
then the approximant will always have a defect at zd, regardless of the degree of [L/P]. 
However, if the pole in f(z) is pole of finite order, it may be possible to remove the defect 
in the approximant by constructing [L/P] with L and P large enough. There is no 
guarantee that such an approximant is possible to construct, but taking higher degree 
approximants frequently removes defects when applied to perturbation theory. 

A defect may also occur in an approximant when f(z) has no poles. In this case, the defect 
could occur due to computational errors which can be significant as discussed in the last 
section. To handle defects of this nature, there are two basic approaches. One approach 
is to use a more accurate method to compute the solution of Equation (3-6). This may 
entail finding the coefficients in the original power series to greater accuracy, using a 
machine capable of higher accuracy or resorting to rational approximants of the 
coefficients of the original power series and then solving Equation (3-6) exactly. The 
second method would be to try calculating the approximant to lower degree in the hopes 
that the C matrix condition number will be enough smaller to prevent the resulting 
approximant from having a defect. 

Another cause for a defect to occur is trying to approximate a complicated function by too 
low an order Fade approximant. If this is the case, merely finding a higher order 
approximant will remove the defect with little effort. Coupling this comment with the last 
point in the previous paragraph, it seems that trying various degree approximants (higher 
and lower) until finding one that is defect-free is a reasonable approach; however, 
checking the condition number of C is an intelligent starting point to determine if the 
difficulty is computational. 

More information on making the calculation of Fade approximants numerically stable can 
be found in Reference [13]. 



I 




APPLYING PADE APPROXIMANTS TO PERTURBATION RESULTS 
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Having briefly explored one type of perturbation theory application, there are some 
obvious opportunities to apply Pade approximants to good ends. 

One of the most direct applications for Pade approximants in vibrations is to predict the 
frequency response of a system in terms of the non-linearity parameter s. There are many 
different methods for estimating the natural frequency of slightly non-linear systems (e « 
1), but many are inaccurate if the system is tuned such that s is slightly larger. Given that, 
it is reasonable to calculate a Pade approximant for the frequency of a non-linear system 
obtained with a perturbation technique in an effort to accurately predict the response 
frequency for larger s. As a starting point, it is assured that a straight power series for ca 
will fail to accurately predict the frequency an oscillatory solution as s increases except in 
very special circumstances. In the general case the frequency will not grow as s" which is 
the asymptotic behavior of the power series for large s. This can be seen by examining the 
power series expansion for the frequency of the van der Pol equation. 

Example 3-4 

Returning to the van der Pol oscillator, Equation (1-1), and using Lindstedt’s method, the 
frequency can be found to high order in s. (This result was derived in Chapter 2.) The 
first few terms are: 



, 1 , 17 , 35 , 678899 , 

= 1-“S + 1 S + - s* +0(s’ ) 



16 3072 884736 5096079360 



Taking the [4/4] approximant: 



1 +- 



[4/4] = 



492389 

2698560 






433361 

34541568* 



1 + 



661049 

^8560 



+ 



1285087 

57569280* 



Graphing both the power series for to (dashed) and [4/4] (solid): 



2 
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Figure (3-3) Power series and Fade approximant to 
the frequency of a van der Pol limit cycle. 



From experimental data, the van der Pol oscillator is in fact a softening system, meaning 
the frequency of the response decreases with increased non-linear effects. Also, just from 
examining the graph, the behavior of the [4/4] approximant seems to continue the pattern 
of the power series expansion, before s became large enough to make © start to grow 
quickly. Both point to the fact that the Fade approximant does a good job of extending 
the approximation for frequency response as 8 increases. 

Large e Concerns 

Another area where Fade approximants are very useful in vibration theory is in application 
to relaxation oscillators. Relaxation oscillators have numerous applications in biological 
systems including how the heart beats, how fish swim (modeling the synapses in the fish’s 
central nervous system), and have possible application as control circuits for robotics. The 
relaxation oscillator is modeled by a non-linear differential equation where the non- 
linearity is allowed to become large compared to the linear components in the system, 
meaning in the neighborhood of 10 or 100 times larger. Consider the relaxation limit for 
the van der Pol oscillator. Equation (1-1), which behaves as a relaxation for large e. As 
seen in Example 3-4, the Fade approximant for the frequency more reasonably predicts 
behavior in the relaxation limit. 

Application to Coupled van der Pol Oscillators 

The main reason for presenting this information on Fade approximants is in applying the 
theory to the stability transition curves that will be derived in Chapter 5. For the 
remainder of this thesis, some basic familiarity with Fade approximants will be assumed as 
they will be used to estimate the frequency of the limit cycle of Equation (1-1), the Zero 
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Mean Damping Surface described in Chapter 4, and the Stability Transition Curves found 
in Chapter 5. 



CHAPTER 4; COUPLED OSCILLATOR ANALYSIS 



PROPOSED MODEL FOR COUPLED OSCILLATORS 

As discussed in the introduction, the pair of van der Pol oscillators (x and y) will be 
modeled by two van der Pol equations coupled with a linear spring and a linear damper. 
The coupling will be “weak” in the respect that the spring and damper will both be order 
s; in particular the stiffhess constant will be sA and the damping constant will be eB where 
8 is the same non-linear parameter as in the original oscillator. The coupled equations take 
the following form [18]: 

X - s(l - )x + X = eA(y - x) + sB(y - x) 
y-e(l-y^)y + y = sA(x-y) + 8B(x-y) 

where dots represent derivatives with respect to time x. 

The object is to determine the existence and stability of solutions to Equation (4-1). In 
particular, it is desired to find what values for sA and eB result in orbitally stable 
solutions. 

In this context, orbitally stable is most easily described in the phase plane. Consider the 
phase plane of the original van der Pol equation plotted in Figure (1-1); the limit cycle 
presented is orbitally stable. In mathematical parlance, for any small perturbation of the 
system away from the described limit cycle, there exists a finite distance 6 such that the 
maximum distance between any point on the disturbed path and the original limit cycle is 
less than or equal to 6 [14]. Understanding the concept of orbital stability is essential to 
comprehending the behavior of coupled oscillators. 

By inspection, an exact in-phase motion of the coupled oscillators is possible with 
x(x)=y(x). In that case, the coupling terms identically vanish and x(x) and y(x) will have 
to satisfy the original van der Pol equation as found in Chapter 1 with x(x) = y(x) = u(x). 
In this situation, the behavior of each of the coupled oscillators is exactly the same as the 
original solution, namely each oscillator would exhibit an orbitally stable limit cycle in its 
own two dimensional phase plane, as seen in Figure (1-1). 
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In the same vein, it is possible to take advantage of the known properties of u(x) by 
considering x(x) and y(x) as functions which vary from u(x) by small amounts. Treating x 
and y as variations away from the solution of the van der Pol equation clarifies the 
meaning of “stability” as applied to the coupled equations of Equation (4-1). While 
studying the variations away from u(x), determining the stability of x and y is equivalent to 
studying the origin in the variational space. This will be addressed in more detail after 
deriving the variational equations. 



DERIVATION OF THE VARIATIONAL EQUATIONS 

Treating x(x) and y(x) as small variations from u(x), the stability of the solutions of 
Equation (4-1) can be determined by characterizing the behavior of the variations. 
Consider both functions as small variations from the limit cycle u(x) [28]: 



^ = x-u => x = ^ + u 

r| = y-u => y = ri + u 



(4-2) 



where ^ is the variation of x(x) from u(x) and r\ is the variation of y(x) from u(x). 
Replacing x and y in Equation (4-1) and rearranging terms yields: 

^ - e(l - - 2u^ - u^)^ + (l + 2 euu)^ + (ii - 8(l - u^)u + u j = 

eA(ti-^) + £B(ti-^) 

(4-3) 

f| - e(i - - 2ut| - u^ )ift + (l + 2Euu)r| + |ii - e(i - u^ )ii + u| = 

£A(^-ti)+eb(4-ti) 



where the terms in braces are the original van der Pol equation and are identically equal to 
zero, resulting in: 



^ - e(i - - 2u^ - u^)^ +(l + 2 euu)^ = eA{ti - ^) + eB(i) - i) 

f) - e(i - T|^ - 2uq - u ^ )tj + (l + 2 euu)t| = eA(^ - ^) + sb(^ - q) 



(4-4). 



In the interest of studying only small variations from the limit cycle, non-linear terms in ^ 
and q are neglected: 



I - e(i - +(l + 2 euu)^ = EA(r|-^) + EB(fi-^) 

ii - e(i - u^)fi +(l + 2 euu)ti = eA(^ - n) +eb(^ - f|) 



(4-5). 
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In order to simplify the coupling terms appearing on the right hand sides of Equation (4- 
5), two new equations can be formed by adding and subtracting the two Equations (4-5) 
and letting: 



w = ^ + q 
v = ^-q 



(4-6). 



So Equations (4-5) become: 

w-E(l-u^)w-i-(l-f-2EUu)w = 0 

( 4 - 7 )' 

V - e(i - 2B - u^ )v + (l + 2 euu + 2eA)v = 0 

By inspecting Equations (4-7) it can be seen that both equations are exact differentials 
which can be integrated with respect to time to yield: 

w-e(1-u^)w-i-w = k, 

/ ^ (4-8). 

V - e(1- 2B - u^)v + (l + 2eA)v = kj 



While studying stability, the solutions for v(t) and w(t) can be recast to further simplify 
Equation (4-6). The particular solutions to Equations (4-8) are constants, in particular: 

Wp =k, 

k, (4-9). 

“ 1 + 2eA 

These solutions do not vary with time and can be removed by a simple coordinate 
transform in the variational space: w-^ w-Wp and v-»v-Vp. Such a transformation will not 
alter the stability of the solutions, since stability is controlled by the exponential behavior 
of the functions. Then the varational equations become: 



w - e(i - u^ )w + w = 0 
v-e(1-2B-u^)v + (1 + 2eA)v = 0 



(4-10). 
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In order to study the stability of the coupled van der Pol oscillators, both of Equations (4- 
10) must be explored and characterized. In order to discover the nature of the first 
equation in Equation (4-10), consider a small variation from a single van der Pol solution, 
u(t); 



z(t) = u(t) + ^(x) (4-11) 

where u(t) is a solution of the van der Pol equation. Equation (1-1). Substituting z into 
Equation (1-1) results in: 

u + ^ + s|l-(u + ^)^|(u + ^) + u + ^ = 0 (4-12) 

which when expanded and eliminating the higher order terms in ^ results in: 

{ii + e(1 - u^ )u + u} + ^ + e(l - u^ )^ + ^ = 0 (4-13) 

where the term in braces is again the van der Pol equation and is identically zero, resulting 
in: 



^+s(l-u^)^ + ^ = 0 (4-14). 

Equation (4-14) is the behavior of a small variation, away from the orbitally stable van 
der Pol limit cycle, therefore ^ must be orbitally stable. Comparing Equation (4-14) to the 
first of Equation (4-10), w and ^ have identical solutions, therefore w(x) must be orbitally 
stable. 

Understanding the stability of the first Equation (4-10) leaves only the behavior of the 
second variational equation: 

v-e(1-2B-u^)v + (1 + 2sA)v = 0 (4-15) 

to be characterized in order to understand the stability of Equations (4- 1 ). Equation (4- 
1 5) will be referred to as the variational equation for the remainder of this thesis. 

SIGNIFICANCE OF THE ZERO MEAN DAMPING SURFACE 

Studying the stability of the variational equation requires characterizing the solutions of 
Equation (4-15) in the three dimensional parameter space { e , eA, eB }. Ultimately, it 
should be possible to find the surface in { e , eA, eB } space which corresponds to the 
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transition from stable to unstable oscillations of the coupled equations. Such an analysis is 
feasible, however it would be formidable to complete using analytical methods. Such an 
analysis using numerical techniques is underway by Reinhall and Storti[l9,20]. 

In order to reduce the complexity of the problem, the parameter space may be reduced to 
two dimensions by fixing or specifying any of the three parameters. Such an analysis for 
the case of B=0, i.e. no damping between the coupled oscillators, is presented by Storti 
and Reinhall [26]. 

Another obvious two dimensional subspace for studying the variational equation is the 
Zero Mean Damping (ZMD) surface of the parameter space. The ZMD surface is the 
locus of values for the damping parameter sB, as a function of s, resulting in zero time 
average damping in the variational equation. The ZMD surface is a logical place to study 
the behavior Equation (4-15) because it marks the boundary between the region of 
negative mean damping, where all solutions to Equation (4-15) will grow over time due to 
the energy added by negative damping, from the region where damping is positive on 
average and stable solutions may be possible. By examining Equation (4-9), the parameter 
sA does not appear in the damping term, so the ZMD surface will be a function of e and 
eB only. The ZMD surface is found by requiring the time average of the coefficient of v 
to vanish in Equation (4-15); 



where the angled brackets indicate time averaging. Clearly Bzmd vvUl be a function of e 
since u(t), found in Chapter 1, is a power series in e. 

The balance of this thesis will be dedicated to characterizing the behavior of the variational 
equation. Equation (4-15), along the ZMD surface. 

CALCULATING THE ZMD SURFACE 

Performing the time averaging using the solution for u(t) found in Chapter 1 and is a 
tedious, but straight-forward process. Since Bzmd is not a function of time. Equation (4- 
16) can be rewritten as: 




(4-16) 




(4-17). 



Using the solution u(t) generated in Chapter 1, the ZMD surface can be calculated to the 
same order in e as u(t). The first few terms of Bzmd are: 
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B 



ZNID 



_ 1 il 48598*^ 12921629s* 

“ ~ 2 ~ 32 9216 10616832 ~ 152882380800 



+ 0(8'”) (4-18). 



A Fade approximant to the ZMD surface through 0(8^®)is given by: 

-0.5 -0.185e2 - 0.0669e'* - 0.01 lOe^ - 0.00132e® - 0.0000589e*0 



Bzmd = 



1. + 0.307e2 + 0.1 17e4 + 0.0163 e6 + 0.00202 e8 + 0.0000654€iO (4-19) 



where the fractions which actually appear in the Fade approximant have been replaced 
with 3 decimals of precision. Figure (4-1) is a plot of the Fade approximant [10/10] of 
Bzmd calculated using exact coefficients. 




Figure (4-1) Fade approximant to the 
Bzmd vs. 8. 
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Figure (4-2) Fade approximant to the 
ZMD surface in 3-D parameter space. 



Figure (4-2) shows no additional information from Figure (4-1), but gives a physical feel 
for the parameter space. The stability transition curves obtained in Chapter 5 should be 
plotted along this surface in 3-space. 

FLOQUET TFEORY AND STABILITY TRANSITIONS CURVES 

Equation (4-15), the variational equation, is a Hill’s equation, meaning it is a linear 
differential equation with periodic coefficients. This has two important consequences in 
relation to examining the stability of the solutions of Equations (4-1) along the ZMD 
surface. 

Since Equation (4-15) is a Hill’s equation whose periodic coefficients have period K, all 
stable, periodic solutions will have period n or 27t. Recalling that u(t) is periodic, with 
period 2ti, u^(t) must have period 7t; then the claim follows directly from Floquet’s 
Theorem [16]. 
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The second result of the variational equation being a Hill’s equation, is that the transition 
from stable to unstable behavior in the parameter space will always be through a stable, 
periodic solution. In combination with the first result, that stable, periodic solution must 
be of period n or lit. To prove this claim, consider Floquet’s Theorem[14]: 

The regular system i = P(t)x, where P is an n x n matrix function with 
minimal period T, has at least one non-trivial solution x = x(t) such that 
x(t + T) = M.x(t) where p is a constant. 

The variational equation can be written in the form required by Floquet’s Theorem by 
selecting; 




and identifying; 



(4-20) 



0 1 
-(l + 2eA) 8(1-2 B^-u^) 



(4-21) 



which has period k due to the presence of u^(t). Therefore Equation (4-15) meets the 
conditions of Floquet’s Theorem. Now consider the possible values for p. If p>l the 
solution x(t) becomes unbounded over time since the solution is multiplied by p for each 
interval of time k; similarly if p<-l the solution will also be unbounded over time. In both 
of these cases the solutions to the variational equation are unstable. If -l<p<l, the 
solution is stable and decays to the origin over time. That leaves two possible values for 
p; p=l or p=-l. 



If p=l, the solution is periodic with period n, the same period as P(t). If p=-l, the 
solution is also periodic but has period 2it. Identifying where the characteristic number, p, 
takes on the values of 1 and -1 marks the transition between stable and unstable solutions 
of the variational equation. 

The fact that the transition between stable and unstable behavior is always through a 
stable, periodic solution with period n or 27c focuses the study of this thesis. Locating the 
periodic solutions of Equation (4-15) on the ZMD surface identifies the stability transition 
curves, where the stability transition curves are the curves which divide stable and unstable 
behavior of Equations (4-1) in the parameter space along the ZMD surface. Hence, the 



28 

next task is to identify where periodic solutions (of period k or 2 ti) of Equation (4-15) 
exist on the ZMD surface. 



CHAPTER 5: PERIODIC SOLUTIONS OF THE VARIATIONAL EQUATION 



As discussed at the end of Chapter 4, identifying the periodic solutions of the variational 
equation along the Zero Mean Damping (ZMD) surface is equivalent to finding the 
stability transition curves for the coupled oscillators modeled by Equations (4-1). The 
next task is to find those periodic solutions. 

SOLUTION METHOD 

The periodic solutions of the variational equation will be identified using Lindstedt’s 
method [14], which was also used in Chapter 1 to solve the van der Pol equation. The 
variational equation along the zero mean damping surface is restated here for easy 
referencing: 

v-8(1-2B-u^)v + (1 + 2sA)v = 0 (5-1). 

Some characteristics of the Equation (5-1) are worth mentioning. The most significant 
feature is that unlike the van der Pol equation itself, the variational equation is a linear 
ordinary differential equation. Since u(x) can be determined to a suitably large order in s 
as outlined in Chapter 1, and Bzmd is the expansion for the ZMD surface, all the 
coefficients in the variational equation are known except sA. As a consequence of being a 
linear equation, solutions to Equation (5-1) will also obey the principle of superposition. 

Even though Equation (5-1) is a linear differential equation, Lindstedt’s method is used to 
find v(t) because of the presence of u(x) which is a power series. In keeping with 
Lindstedt’s method, the first step is to rescale time such that x = (o(s) t. Substituting this 
rescaling into Equation (5-1) and changing the differentiation to the new time scale; 

v"-s(i-2B2md -u^)o)v' + (a^(l + 2EA)v = 0 (5-2) 

where primes represent differentiation with respect to the rescaled time t. Notice that here 
the function ©(s) was already determined in solving Equation (1-1) and is given by 
Equation (1-19). The fi'equency expansion must be the same as for the original van der 
Pol solution because the solutions to Equation (5-1) must have the same period or twice 
the period of its coefficients; this was discussed at the end of Chapter 4. Plugging in 
power series expansions for o, v, Bzmd, u and sA according to: 
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sA = Ao + J]Ai8' 

i=l 

n,2 

®ZN1D ~ 

i=2 

u(t) = Uo(t) + 2ui(t)e' (5-3) 

i=l 

v(t) = Vo(t) + Xvi(t)s' 

i=l 

n,2 

(0 = 1 + 

i=2 

and collecting powers of e gives a series of linear, ordinary differential equations, the first 
three of which are presented here: 

v" + (l + 2Ao)Vo = 0 

v"+(l + 2Ao)v, =-2A,Vo+v;-2BoV;-uJv; (5-4). 

v" + (1 + 2 Ao)v2 = - 2 A 2 V 0 - 2A,v, - 2uoU,v; + v,' - 2Bov| - uj vj - 2o ju" 

Notice that the left hand side of each of Equations (5-4) is of the same form, namely: 

LHS = v"+(l-^2Ao)Vi (5-5). 

Insisting that the solution of Equation (5-1) be periodic with period 7t or 2k, as discussed 
above, the coefficient for Vj in Equation (5-5) must be an integer squared. This follows 
directly from the Floquet theory discussed at the end of Chapter 4; the solution to 
Equation (5-5) is Cos(t^l - 2A(, ) which can produce v(t) with period K or 2n only if 

VI- 2 A 0 is an integer; therefore: 



-1 

l + 2Ao = k^ or Aq = — - — where k e {l,2,3,...} (5-6). 

This gives rise to an infinite set of starting values for Ao, and possible stability transition 
curves out of each starting value. The first few Ao’s are given here: 




(5-7). 



DEALING WITH Ao=-l/2 
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When Ao = -1/2, the variational equation changes form significantly, the zeroeth order 
equation of Equations (5-4) becomes: 

v" = 0 (5-8) 

which has solution; 

Vq = +Co,t (5-9) 

However, since the goal is to find periodic, stable solutions to the variational equation, Coi 
= 0 so Vo = Co is the solution. Plugging in this Vo and examining the first order equation 
yields: 



v" = -2A,Co (5-10). 

In order for vi to be bounded, either Co = 0 or Aj = 0. The former yields a trivial order 
zero solution, but the latter is a valid method of solving Equation (5-8) wdth a bounded 
solution. By similar arguments used about Equation (5-9) it can be seen that Vi = Ci. 
Following this pattern, the same equation arises at each order of s, so that the stability 
transition curve from Ao=-l/2 is in fact just sA= -1/2. 

SOLVING FOR THE STABILITY TRANSITION CURVES 

For the remainder of the possible values for Ao listed in Equation (5-7), the solution 
procedure is a straight-forward application of Lindstedt’s method. Whereas the solution 
of the van der Pol equation in Chapter 1 resulted in fixing the frequency of the response as 
a function of e, solving the variational equation for v(t) will give the coefficients for the 
sA expansion. 

Unfortunately, unlike the arguments made in Chapter 1 for the van der Pol equation, there 
are no obvious symmetry arguments to use in the solution of the variational equation. 
While the solution of the variational equation shows interesting patterns for some values 
of Ao, the overall pattern is not easily discernible. This results in a very costly 
computational procedure for determining the stability transition curves. 

The stability transition curves calculated below and Tabulated in Appendix E were 
generated using the Mathematical code found in Appendix B. One feature of the code is 
worth mentioning here. In an effort to make solving the second order differential 
equations more efficient, a specific solving routine called MySolve[ _ ] was written. This 
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routine takes advantage of the fact that at every order of s, the differential equation 
encountered always takes the form: 

v"+k^V; = f(sin(mt),Cos(mt)) wherem= 1, 2, 3,...,2i+l (5-11) 

and allows the code to produce solutions without calling a general differential equation 
solver. The specialized solver significantly increases the computational efficiency. 

CONDITIONS FOR SUPPRESSING SECULARITY 

Following Lindstedt’s method and substituting the power series of Equation (5-3) into 
Equation (5-2) and collecting terms of the same order in s, the first few equations can be 
obtained: 

k^vfO, t] +v®’2>[0, t] = 0 

k^v[l, t] t] = -2a[l] v[0, t] +2v®’”[0, t] -4Cos[t]^ v®’”[0, t] 

k^\[2, t] + v®’ 2)[2, t] = -2 a(2] v(0, t] - 2 a[l] v[l, t] - 3 Cos[t] Sin[t] v®>‘>[0, t] + 

Cos[t] Sin[3 1] t] +2 t] -4Cos[t]^ t] + - v®>^>[0, t] 

(5-12) 

where k^ = 1 + 2 Ao as defined in Equation (5-6). The notation used in Equations (5-12) 
is consistent with the Mathematica™ code found in Appendix B. In this notation, v[i,t] is 
the same as Vj(t), a[i] is Ai, and the numbers inside parenthesis represent differentiation 
with respect to the variable t. 

Noting that the left hand side of each of Equations (5-12) are identical, it can be seen that 
the homogeneous solution for each order of E will be similar: 

vh[i, t] = c[i] Cosfkt] + s[i] Sin[kt] (5-13) 

where c[i] and s[i] will be used consistently to represent the coefficients of the 
homogeneous solutions of the f** order equation. The particular solutions will not be as 
consistent, but will be of the form. 

2i+l 

Vp[i, t] = ^ 1 c[i, m] Cos[m kt] + s[i, m] Sin[m k t] } 

m=o (5_14). 



The c[i,m] and s[i,m] are tabulated in Appendix D. 
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In the course of solving Equations (5-1 1), two unique solutions for suppressing secularity 
will arise for each Ao [20], The outcome of these two possible solutions for v[t] and eA is 
that for each starting value of Ao, other than -1/2, there will be two curves along the ZMD 
surface which have valid, periodic solutions, with period 7i or 2 ti. These curves intersect 
at eA = Ao and resemble a horn with its point downward and will sometimes be referred 
to as the “stability horn” arising from each Ao. Extending this result into the three 
dimensional parameter space of Figure (4-2), the two dimensional surface where periodic 
solutions to Equation (5-2) exist will resemble cones will their points coming out of each 
Ao. The regions within the cones are unstable, that is the solutions to the variational 
equation will grow over time, making the origin of the variational space an unstable 
equilibrium point. The regions outside of the cones and above the ZMD surface will be 
asymptotically stable; all solutions will decay towards the origin of the variational space 
over time. The actual surfaces of the cones are configurations which yield orbitally stable 
solutions. 

One of the main difficulties in solving Equations (5-12) will be to identify this split in the 
stability conditions. The conditions are different for each Ao and will be considered one at 
a time. 

Case 1; Ao= 0 

For the case of Ao = 0, Equation (5-6) yields k=l. Solving Equations (5-12) and 
suppressing the secular terms appearing on the right hand side of the first order term 
gives; 



-c[0]-2a[l]s[0]=0 
-2a[l]c[0] -s[0] =0 

which has two possible solutions for a[l] and s[0]: 

({s[0] -» -c[0], all] 1), {s[0] -» c[0], a[l] -y)) 

The first solution will yield a curve to the right of the second when viewed on the ZMD 
surface. From this point on, the solutions for v[t] and eA are unique and the stability 
transition curves are given by; 



(5-15) 

(5-16). 



^ow = 


e 


"1; 

1 


3e5 


247^ 


11657e’ 


1224811 


32 


384 1024 


442368 


* 21233664 


1 

2548039680 


eAup = 






3e5 




1 1657 e’ 


122481 le« 


32 ~ 


384 


1024 


1 

442368 


1 
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1 
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+ 0(e^) 
0(e^) 



(5-17). 
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The power series given by Equation (5-17) were calculated out to 0(s^^) and converted 
into Fade approximants which are plotted in Figure (5-1) with the power series. 




Figure (5-1) Power series and Fade approximants 
[16/16] to transition curves out of Ao=0. 



Case 2: Ao= 3/2 

For Ao = 3/2, Equation (5-6) gives k=2. Suppressing the secular terms in the first order 
equation of Equation (5-12) yields; 



- 2 a[l] s[0] =0 

-2a[l]c[0] = 0 (5-18) 

which requires a[l]=0 for a non-trivial solution. Suppressing secularity in the second 
order equation: 



2 s[0] 
3 

c[0] 

3 



-2 a[2] s[0] = 0 
-2 a[2] c[0] =0 



which has three possible solutions; 

{{a[2] -i s[0] -» o}, (a[2] -» j, c[0] o}, {s[0] 0, c[0] -» 0}} 



(5-19) 



(5-20). 
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The last solution is the trivial solution, but the first two again provide the split for the 
lower and upper stability curves. From this point, the solutions are unique and the stability 
transition curves are given by; 

3 109^4 3419 e6 297305 c* _ ,o 

tA(o„ — + + O (e ) 

2 6 3456 1990656 573308928 

3 53^ 6983 740213 c* _ 

^23 3456 1990656 2866544640 (5-21). 

The Fade approximants for the power series of Equations (5-21) calculated to 0(e^^) are 
plotted in Figure (5-2) with the power series. 




Figure (5-2) Power series and Fade approximants 
[12/12] to transition curves out of Ao=3/2. 



Case 3: Ao = 4 

For Ao = 4, Equation (5-6) gives k=3. Suppressing the secular terms in the first order 
equation of Equation (5-12) yields; 

- 2 all] s[0] =0 

-2a[l]c[0] = 0 (5-22) 



which requires a[l]=0 for a non-trivial solution. Suppressing secularity in the second 
order equation; 
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9s[0] 

16 

9c[0] 

16 



-2a(2] s[0] = 0 
■ -2a[2] c[0] =0 



(5-23) 



which requires a[2] — 9/32 for a non-trivial solution. Suppressing secularity in the third 
order equation: 



5 c[0] 

a[3] s[0] = 0 
5 sfO] 

-2 a[3] c[0]+— ^=0 
32 

which has two possible solutions for a[3] and s[0]: 
{{s[0] -c[0], a[3] (s[0] c[0], a[3] 



(5-24) 



(5-25). 



These two solutions provide the split for the lower and upper stability curves. From this 
point, the solutions are unique and the stability transition curves are given by: 



. ^ 9 5 2877 

— 4 + + 

32 64 40960 

2537 547076^ 11 663959 



737280 15728640 



14155776000 
2537 



9 ^ 5 ^ 2877 

cAm = 4 1 + 

^ 32 64 40960 737280 



17260922723 ^ 
12683575296000 



+ O(^) 



1 7260922723 „ 

h O (t ) 

12683575296000 

547076^ 11663959 6^ 

15728640 ~ 14155776000 



(5-26). 



The Fade approximants for the power series of Equations (5-26) calculated to 0(e^^) are 
plotted in Figure (5-3) with the power series. 
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Figure (5-3) Power series and Fade approximants 
[16/16] to transition curves out of Ao=4. 



CHAPTER 6: SUMMARY AND FUTURE WORK 



In Chapter 5, the first three pairs of stability transition curves along the Zero Mean 
Damping (ZMD) surface were determined. Combining Figures (5-1,2 and 3) yields the 
stability regions for the smaller values of sA. In Figure (6-1), U marks regions that are 
unstable and S marks regions where orbitally stable solutions to the variational equation. 
Equation (4-15), exist. The stability transition curves are approximated by Fade 
approximants. 




Figure (6-1) Stability regions for coupled 
VDP oscillators on the ZMD surface. 

The three dimensional version of the same graph is given in Figure (6-2) where the three 
dimensions are the three parameters of Equation (4-1), (e, eA, eB}. This is a graph of the 
stability transition curves found in Chapter 5 plotted on the ZMD surface of Figure (4-2). 

SUMMARY 

Using the techniques outlined in Chapter 5 and employing a symbolic mathematics 
package (Mathematica™) the first three pairs of stability transition curves for a pair of 
weakly coupled van der Pol oscillators were determined. All calculations were carried out 
exactly, that is exact fractions were used in all computations with no numerical roundoff. 
Using this analytical approach, the six stability transition curves seen in Figures (6-2) and 
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(6-2) were calculated to 0(s^^). The power series coefficients for the transition curves are 
tabulated in Appendix E. 

FUTURE WORK 

There are three main areas for future work on this project. 

The first is to extend the range of sA covered in this thesis, namely find the stability 
transition curves coming from the rest of the Ao values given in Equation (5-7). Work on 
this area, using the code written for this thesis, is already in progress. The main difficulty 
in finding the stability transition curves coming from the next few values of Ao is showing 
that the curves are unique. Present conjecture is that there are more conditions needed to 
fully specify the upper curves than used in Chapter 5. To find what those conditions might 
be, research is underway to use Hill’s determinants to find the periodic solutions rather 
than Lindstedt’s method. 

Another area for future work was hinted at in Chapter 4 when the ZMD surface was 
introduced. The full challenge for this type of stability analysis would be to find the 
stability transition surfaces in the three dimensional parameter space. Currently, numerical 
calculations of this project are under way by Storti and Reinhall (REFERENCE>»). The 
logical extension of this project would be to find the transition surfaces analytically to 
avoid numerical roundoff errors. 

The final, and ultimate continuation of this work is to physically realize the coupled 
relaxation oscillators modeled by Equation (4-1). Once these oscillators are accurately 
predicted by the mathematical models, they could be used to generate control signals for 
ambulatory robots, mechanical fish or similar mechanisms. The goal of this thesis was to 
further understanding toward reaching this goal. 
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Figure (6-2) Three dimensional plot of the 
stability transition curves along the ZMD surface. 
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APPENDIX A: MATHEMATICAtm CODE FOR VDP EQUATION 



This first block of code produces the first three terms of the of the VDP equation. 
Comments are included to help clarify the procedure. A second set of code is included 
which will start with a certain number of terms already in the solution and calculate two 
more terms. That code can be used to extend the solution as far as the user wants or until 
the computer runs out of memory. 

STARTER PROGRAM 

$Line = cd; 

SetDirectory[ ”c:\\vi%>\\catplex\\oda1:a"] 

(★ pinult[f,g,x^n) multiplies power series f by power series g 
where x is the variable, and ke^s only throu^ order n *) 

n n-i 

F»™ult[f_, g_, x_J5ynfcol, n_JEnteger] := Coefficient(f , x, i] Coefficient [g, x, j] 

<* getX[n) creates a power series in e of order n with argument x[i,t] *) 

n 

getX[n_] := E:qDand[^e^x[i, t] ] 



{* xeven[] and xodd[] assign x[i,t] to be a complex cosine or sign *) 



:= ^ A[i, j] (Z^ + Z'^) 



A >2 




>1 
A j»2 



(* getW[n) creates the frequency power series *) 



n 




ic2 
A 1«2 



(★ noslope [] enforces the boundary condition in Equation (1-1) ★) 

noslcpe[i_] := Module [ {tenp } , teftp= Coefficient[X(t] , e, i] ; 

Solve[ (inyD[tenp] / . Z-> 1) == 0, A[i, 1] ] ] 

(* myD[] and myD2[] are specialized derivatives *) 

:= f /. {Z®-- * laZ®, Z-. IZ) 

inyD2[fj := myD[iryD[f] ] 



ClearAll[x, w, X, W, tenp^ Ihs^ rhs. A, B, i] ; (★ start with clean slate *) 
n= 3; (★ order of solution ★) 

X[t] = getX[n] ; (* build power series in e ★) 

W= getW[n] ; 

(★ build the right hand sides of Eq 1-5 *) 

rhs= pmilt[e, Finult[W, prrult[dtX[t] , 1 - piult[X[t] , X[t] , e, n - 1] , e, n- 1] , 
e, n-1] , e, n] ; 

lhs= Fniult[Finult[W, W, e, n] -1, d{t, 2 }X[t] , e, n] ; 

(* assign x[i,t] fom to match Eq 1-15 ★) 

Table[x[i, t] = xeven[i] , {i, 0, n, 2}] ; 

Table[x[i, t] = xodd[i] , {i, 1, n, 2)] ; 

(★ back to building right hand sides *) 

mid= rhs- Ihs / . Table[x^^'^^ [i, t] ni^[x[i, t]] , {i, 0, n}] ; 
inid= mid/ . Table[x^^'^^ [i, t] -> in^)2[x[i, t] ] , {i, 0, n} ] ; 

X[t] = X[t] ; 

(★ pull out each right hand side for each order of e ★) 

tenp= Table[Coefficient[mid, e, i] , {i, 1, n}] ; 



(★ free up some memory *) 

ClearAll[lhs, rhs, mid, x] ; 
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(* Solving the 1st order equation *) 

i= 1; 

ric^it= Collect[Expand[ j ^ Table[Z^^*^, {jr 0/ i)]] / 

sec= E3q>and[Coefficient[ricJit, Z] ] ; (* Identifies the secular terms *) 

soIKule= Solve [ sec == 0, A[0, 1]] ; (★ Suppresses secularity ★) 

A[0, 1] = A[0, 1] /. solRule|[3l ; 

ricjit = rigjit / . solRulei3] ; 

_ Coefficient[rig^t, , 

T^le[A[i, 2 j + 1] = — — jyj , {j, (* Enforces Eq 1-15 *) 

A[i, 1] = A[i, 1] /. noslcpe[i] 111 ; (★ Enforces boundary condition *) 

tenpli] :: 0; 



(* Solving the 2nd order equation *) 



i= 2; 

ri^t= CoUect[E3qpand[tenpIiI] , Table[z^^^^, {j, 0, i)]] ; 

soU^e= Solve [ Expand [ Coefficient [ricjit, Z] ] == 0, w[i]] ; (* Id* s the secular terms ★) 
w[i] = w[i] /. solE?ule|ll ; 
rigjit = rigjit / . soU^elll ; 



Table[A[i, 2 j + 1] 



Coefficient [ricjit, 
1- (2j+l)2 



, {j, 1, i)] ; 



(★ Enforces Eq 1-15 *) 



tenplil = 0; 



(★ Solving the 3rd order equation ★) 
i=3; 

ricjit = Collect[ Expand[ 



] , Table[z2^*\ {j,0, i)]]; 



solI^e= Solve[Expand[Coefficient[rigfit, Z] ] == 0, A[i- 1, 1]] ; (* Id* s the secular temns * 
A[i- 1, 1] =A[i- 1, 1] /. soU^eHl; 
rig^it = rigjit / . solRulelll ; 

r . « . ^ Coefficient [ri^t, i ^ -i -ic 

Table[A[i, 2 j + 1] = , {j, 1, x} ; (* Enforces Eq 1-15 *) 

1- (2 j+ 1)2 

A[i, 1] = A[i, 1] / . noslope[i] [11 ; 
tenplil = 0; 



(★ Enforces boundary condition *) 
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(* Saving the data — *) 

DeleteFile [ "af lie . txt" ] 

Save [ "af ile . txt" , A] 

DeleteFile [ " wf lie . txt" ] 

Save ["wfile.txt", w] 

DeleteFile [ "nfile . txt"] 

Save ["nfile.txt", n] 

(★ This section puts the coefficients back in a 
form to use with sines and cosines . . . * ) 

Table[Table[B[i, j] = 2A[i, j] , (j, 1, 2i+ 1, 2}] , {i, 0, n, 2}] ; 
T^le[Table[B[i, j] = -2A[i, j] , {j, 1, 2i + 1, 2}] , {i, 1, n, 2>] ; 
Save ["bfile.txt”, B] 

(* Rebuilding X[t] in terms of sines and cosines ★) 

2U1 

Tahle[x[i, t] = ^ 2}] ; 

A >2 

2U1 

Table[x[i, t] = ^ B[i, j] Sin[jt], {i, 1, n, 2}] ; 

A >2 

X[t] = Sum[x[i, t] e^, {i, 0, n}] ; 

(* Check out the work... *) 

W 

X[t] 



CONTINUATION PROGRAM 



The continuation program uses the same functions defined above, and the code is identical 
down through defining noslope( ] then change the code as follows; 
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ClearAll[x, w, X, W, tenp^ Ihs, rhs, A, n] ; 
n« "nfile.txt"; 
old= n; 
n = old+ 2 ; 

Print[ "Starting, old = ", old, " new = ", n] 
tl = Tiin^sed[ ] ; 

X[t] = getX[n] ; 

W= getW[n] ; 

rhs= l-Fmalt[X[t] , X[t] , e, n-1] ; 
tenp= pinult[W, dtX[t] , e, n- 1] ; 
rhs = pxultlrhs, tenp, e, n- 1] ; 

Clear [tesip] ; 

rhs= rhs/. Table[e^-+ a^, {i, old, n- 1}] ; 

rhs = rhs / . 0; 

rhs = rhs / . a-+ e; 

rhss pinLilt[e, rhs, e, n] ; 

lhs= piniLt[W, W, e, n] ; 

lhs= Ihs- 1; 

lhs= praalt[lhs, d[t,2) X[t] , e, n] ; 
lhs= Ihs/. TabLe[e^-*a^, {i, old+ 1, n}] ; 
lhs= Ihs/ . e-> 0; 

Ihs = Ihs / . a-* e; 

Table[x[i, t] = xeven[i] , {i, 0, n, 2}] ; 

Table[x[i, t] = xodd[i] , {i, 1, n, 2)] ; 

ndd= rhs- Ihs /. Table [i, t] -► itiyD[x[i, t]] , {i, 0, n}] ; 
mid= mid/. Table[x^°'^ [i, t] -►m^D2[x[i, t]] , {i, 0, n)] ; 
Print ["updating X[t] . . ."] 

X[t] = X[t] ; 

Print[ "building tenp. . . "] 

tenp = {mid / . {e^^^ -► 1 , e 0} , mid / . {e^^^ 1 , e -► 0} } ; 

ClearAll [ Ihs , rhs, mid, x] ; 

A« "afile.txt"; 
w « "wfile.txt"; 



48 



i = old+ 1 ; 

rhs = Collect [ Ejqpand[ tenpllj ] , Table[Z^ , { j , 0 , ^) ] ] ' 



solf^e= Solve [E3q>and[ Coefficient [rhs, Z] ] == 0, w[i]] ; 
w[i] = w[i] /. solRulellJ ; 
rhs = rhs / . solRule[l] ; 



Table[A[i, 2 j + 1] 



Coeffici^t[rhs, 
1- (2j+l)2 



{j. If i>] ; 



terpllj = 0; 
i = old* 2 ; 

rhs= Collect[ Expand[ ] , Table[z^^*^, {j, 0, i}]] ; 



solRule= Solve[Expand[Coefficient[zhs, Z]] == 0, A[i- 1, 1]] ; 
A[i- 1, 1] =A[i-l, 1] /. solRulelll; 
xhs = xhs / . solRule|l] ; 

, , r . « . . Coefficient [rhs, . , 

Table[A[i, 2 j + 1] = _ ^2 j + 1 ) 2 ' 1 



A[i, 1) = A[i, 1] / . noslcpe[i] (Ij ; 
tenp|[2] = 0; 

DeleteEi.le[{ "afile.txt", "wfile.txt", "nfile.txt"}] 
Save [ "af ile • txt" , A] 

Save ["wfile.txt", w] 

Save ["nfile.txt", n] 



From here the last few routines in the starter program can again be used to translate the 
solution back into sines and cosines and save the results in that form. 



APPENDIX B: STABILITY TRANSITION CODE 



The first program is “Recover.nb” which converts the coefficients of the complex VDP 
solution into real coefficients for use in the stability routines. It also calculates the ZMD 
surface coefficients. It updates these coefficients in the c:\stabil\data folder as well as the 
frequency coefficients. This code needs to be run only once after generating as many 
terms in the VDP solution as necessary. 



RECOVER.NB 

Directory structure is assumed to be 
c: 

-complex THIS HAS THE VDP SOLUTION STORED IN IT 

- stabil DATA FOR THE STABILITY CURVES 

$Line= <»; 

SetDixectory["c:\\v«:%>"] ; 

n« "ccnplex\\c3data\\nf ile . txt" ; (* nuniser of teams available in VDP solution *) 

Reads in solution to complex solution and converts them to coefficients for real solution 

A« "oatplex\\cxJata\\af ile . txt" ; 

T^le[ Table [B[i, j] = 2A[i, j] , {j, 1, 2i + 1, 2}] , {i, 0, n, 2}] ; 

Table [ Table [B[i, j] = -2A[i, j] , {j, 1, 2i + 1, 2)] , {i, 1, n, 2}] ; 
DeleteFile["ocnplex\\cx^ta\\bfile.txt"} ; 

DeleteFile[ "st^il\\data\\bfile.txt"] ; 

Save ["caiplex\\odata\\bf ile. txt", B] ; 

Save["stabil\\data\\bfile.txt" , B] ; 

The next four lines calculate the coefficients of the ZMD surface 

2U1 1 

T^le[xli] = ^ - B[i, j] (Z5+ Z-3) , {i, 0, n, 2)] ; 

3-1 ^ 

A 

2i+l 1 

Table[x[i] = ^ — IB[i, j] (Z^-Z'^) , {i, 1, n, 2}] ; 

3-1 2 

4 3-2 

1 1 

T^le[b[i] = ^ |Ebq>and[ x[j] x[i- j]] / . Z — * O) , {i, 0, n, 2)] ; 

3-0 2 

b[0] = b[0] + — ; 

2 

Saves the ZMD coefficients, updates the frequency coefficients in the stabil folder and 
updates how many terms are available in the VDP for the stability routines 
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DeleteFile [ ”stabil\\data\\crit . txt" ] 

Save [ *' stabil\\data\\crit . txt " , b] 

DeleteFile [ "st^il\\data\\wf ile . txt"j 

CcpyFile[ "ccnplex\\odata\\wfile . txt " , *'st^il\\data\\ wfile . txt"] 
DeleteFile [ "st^il\\data\\nfile. txt"] 

Save ["st^il\\data\\nfile. txt", n] 

The next program will start the solution for the stability transition curve: 



STARTER.NB 



Assumes the same directory structure as Recover, nb 

$Line= oo 



SetDirectory( "c:\\v^%>\\st^il\\data"] 

pmult[ ] is the same as in the VDP code. getV[ ] is equivalent to getX[ ] from VDP and 
get A[ ] makes the power series for the transition curve. 

n ivi 

pniult[f__, x_Syn4x>l, n__Integar] Coefficient[f , x, i] Coefficient[g, x, j] 

i«0 >0 



OCdET 

getV[order_j := ^ e^v[n, t] 

nsO 

gatA(order_j := ^ e^a[n] 

HbO 



update[ ] updates the RHSs of Equation (5-12) as each v[i,t] is determined. 

v.pdate[tenp_, i_] := Modole[{ttt) , ttt= tenp; 

tttli- 11 = 0; 

Table[ttt|IiI = tttlH / . v[ j, t] -♦ v[ j, t] , { j, 0, i - 1)] ; 

Table[tttlil = tttlil /. t] -»atv[j, t] , {j, 0, i-1}] ; 

Table[tttlil = tttlij /. t] -»a{t, 2 )V[j, t] , {j, 0, i-1}] ; 

ttt] 

Specialized ODE solver discussed in Chapter 5. 
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iTYSolve[i_, ri^t_] := Mcx±ile[{n, tenp) , 
n= 2 (i+ 1) + 1; 

tenp= c[i] Cos[k t] + s[i] Sin[k t] ; 



k^l 



tenp= tenp+ ^ 



Coefficient [rigfit, Cos[mt]] Cos[mt] 
- m2 + k2 



A Coefficient[ri^t, Cos[m t] ] Cos[m t] 

tenp= tenp+ > 

mflEli -m2+k2 

Coefficient [ri^t, Sin[mt]] Sin[mt] 
tesip=tenp+ > 

m=l 



-m2+ k2 



^ Coefficient [ri^t, Sin[mt]] Sin[mt] 
tesip= tenp+ > 

ZllslCfl 

tenp] 



- m2 + k2 



Build the power series for the frequency, W, the VDP solution, U, and the ZMD surface, 
B. 

5 

W= 1+ 2 w[i] ; 

A i»2 
5 

t]; 

i=0 

5 

B= 2 e^b[i] ; 

ilO 
A U2 

ClearAll[n, v, V, A, a, tenp, r±<^t, s, c] 

n= 5; (* how many terms to calculate *) 

V[t] = gstV[n] ; 

A= g^tA[n] ; 

a[0] =0; (* v^ch val\je of a[0] to go after *) 

k = V 1+ 2a[0] (* k of Eq (5-6) *) 



rhs is the right hand side of Equation (5-12) prior to separating the powers of e. 

Once they are separated, each power is stored in temp[i]. 

rhs= pnult[pinL2lt[eW, ^tV[t] , e, n] , 1 - pinLilt[U, U, e, n- 1] -2B, e^ n] - 

2Finult[A- aO, V[t] , e, n] - pirult[pinLilt[W, W, e, n] - ^{t, 2 > V[t] , e, n] ; 

tenp a Table [Coefficient [rhs, e, i] , {i, 1, n}] ; 

Clear [rhs] ; 



Read in the coefficients for VDP, frequency and ZMD surface (generated by recover. nb). 
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u<< ’’newu.txt"; 
w<< "wfile.txt”; 
b<< ”crit.txt”; 

Solve the zeroeth order equation, note the rhs is zero for this equation. Last line is the 
condition of unit displacement for the perturbation. It can be replaced with other 
conditions or left out entirely. 
i= 0; 

v[i, t] = n^solveti, 0] ; 

v[i, t] = v[i, t] /. Solve[<getV[i] / . t-» 0) == 1, c[i]]IU; 



Solve the first order equation. First update the rhs to reflect v[0,t], next find the secular 
term coefficients, si and cl. Next solve to eliminate the secular terms, in this case for a[l] 
and s[0]. Next update the coefficients and the rhs with the solution just found and finally 
call the specialized ODE solver to get v[I,t]. Again the last line imposes unit displacement 
and could be left out or replaced by other conditions. 

Note; For this case, a[0]=0, picking sol[[l]] gives the lower curve while sol[[2]] would 
give the upper curve. 
i= 1; 

tenp= i;^pdate[tenp, i] ; 

ri^t = Ck3llect[Tric^<educ3e[tenp[i] ] , (Sin[ t] , Cos[t] } ] ; 

si = Coefficient [ri^t, Sin[t] ] ; 

cl = Coefficient [ri^t, Cos[t] ] ; 

sol = Solve[{sl== 0, cl== 0} , {a[i] , s[i- 1]}] 

a[i] = a[i] /. solfl] ; 

v[i-l, t] =v[i-l, t] / .sol[l]; 
s[i- 1] =s[i- 1] /. soUll; 
ri^t = ri^t / . sol[l] ; 
v[i, t] = n^sol've[i, ri^t] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] / . t-» 0) == 1, c[i]]|[l]l; 

Solve second order equation, just like above. 

i = 2; 

tenp= vpdate[tenp, i] ; 

ri^t= Collect[Tri g R e duo e [ tenp|il ] , (Sin[t] , Cos[t]}] ; 

si = Coefficient [ri^t, Sin[t] ] ; 

cl = Coefficient [ri^t, Cos[t] ] ; 

sol= Solve[(sl == 0, cl == 0) , <a[i] , s[i- 1] }] 

a[i] =a[i] /. sol[l]; 

s[i- 1] = sti - 1] / . sol[l] ; 

ri^t = ri^t / . sol[l] ; 



v[i, t] = iriysolve[i, ri^t] ; 

v[i, t] = v[i, t] /. Solve[(getV[i] / . t-> 0) == 1, c[i]]|[l]|; 

Same thing for third order. . . 

i= 3; 

tenp= -updatGltesqp , i] ; 

ri^t= tollect[Tricpedbc3e[ , {Sin[t] , Cos[t] }] ; 

si = Coefficient [riqjit, Sin[t] ] ; 

cl = Coefficient [ri^t, Cos[t] ] ; 

sol= Solve[{sl== 0, cl== 0}, {a[i] , s[i- 1]}] 

a[i] = a[i] /. solllj ; 

s[i- 1] = s[i - 1] / . sollll ; 

right = rigjht / . sol|[l]| ; 

v[i^ t] = irysolveti^ ri^t] ; 

v(i, t] = v[i, t] / . Solve[ (getV[i] / . t-» 0) == 1^ c[i] ] 111 ; 



Same thing for fourth order... 

i= 4; 

tenp = \;pdate[tenp^ i] ; 

ri^t= Collect[T 2 i.^teduoe[tenp[il ] , {Sin[t] , Cos[t]>] ; 

si = Coefficient [rights Sin[t] ] ; 

cl = Coefficient [ri^t, Cos[t] ) ; 

sol = Solve[{sl == 0^ cl == 0} , {a[i] , s[i - 1] }] 

a[i] = a[i] /. soli 11 ; 

s[i-l] =s[i-l] /. sollll; 

rigjit = rig^it / . sollll ; 

v[i^ t] = nysolve[i, ric^it] ; 

vti, t] = v[i, t] / . Solve[ (getV[i] / . t-^ 0) == 1, c(i] ] HI ; 

Same thing for fifth order... 

i= 5; 

tenp= i:pdate[tenp^ i] ; 

right = Collect [Tricpedbc»[tenp|[ii ] , {Sin[t] , Cos[t]}] ; 

si = Coefficigit[ri^t, Sin[t] ] ; 

cl = Coefficient [ri^t, Cos[t] ] ; 

sol= Solve[{sl == 0, cl == 0} ^ {a[i] ^ s[i - 1] }] 

a[i] = a[i] /. sollll; 

s[i- 1] = s[i - 1] / . sollll ; 

ri^t = rigjit / . solHl ; 

v[i, t] = irYsolve[i, right] ; 

v[i, t] = v[i, t] /. Solve[(getV[i] / . t-^ 0) == 1, c[i]]lli; 
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Save the results. aOl.txt is the lower curve out of a[0]=0. nOl.txt is how many terms of 
v[i,t] have been solved. 

DeleteFile [ "aOl . txt" ] 

Save [ "aOl . txt " , A] 

DeleteFile [ "vOl . txt” ] 

Save["v01.txt", V] 

DeleteFile [ "nOl . txt" ] 

Save ["nOl. txt", n] 



CONTINUE.NB 

Continue \vill take the results of either starter.nb, or itself and extend the results by four 
terms. The functions down through and including mysolve[ ] from starter.nb need to be 
included as well as MyExpand[ ] which uses Euler’s identity to update the right hand side 
of Equation (5-12) more quickly: 

l^^E: 7 >and[insss_J := Mrxiile[{Z, tenp} , 

tenp= Expand[mess / . (Cos[m_t] -> 1 (Z”*+ Z""*) , Sin[m_t] -* -1 I (Z™- Z"™) , 

Cos[t] ^ - (Z. -) , Sxn[t] - -- I (Z- -)}] ; 

E:q>and[tenp / . — ► Cos[m t] + I Sin[m t] , Cos[2 1] + I Sin[2 t] , 

Z-^ Cos[t] + I Sin[t] }] ] 

The next section of code checks to ensure enough terms are available in the VDP solution 
and then builds the right hand sides of Equation (5-12) much as above. It prints out lines 
to let the user know what order terms it is trying to calculate: 

n« "nfile.txt"; 
m= n; 

Print[m, " terms in ^^P solution available."] 
n<< "nOl.txt"; 
min= n+ 1; 
n= min+ 3; 

Print ["Starting with ", min, ", going to ", n, "."} 

If[n> m, Print["Not encu^ termsl i"] ;Quit[] , Clear[m]} 
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W= 1+ ^ w[i] ; 

i=2 
A U2 

n 

U= ^e^u[i, t] ; 

i-0 

n 

ZM3= 2] e^b[i] ; 
feo 

A is2 

ClearAll[v, V, A, tenp, rhs, s, c] 

V[t] = getV[n] ; 

A= getA[n] ; 

rhs= pinult[eW, dtV[t] , e, n] ; 

tenp= l-Finult[U, U, e, n] -22MD; 

rhs= Fitult[rhs, tenp, e, n, min] ; 

tenp= 2pimalt[A- a[0] , V[t] ^ n, min] ; 

tenp2 = pnraLilt[W, W, e, n] - 1 ; 

tenp2 = pinult[tenp2, d{t, 2 > V[t] , e, n, min] ; 

rhs = rhs - tenp - tenp2 ; 

ClearAll[tenp, tenp2] ; 

tenp= Table [Coefficient [rhs, e, i] , {i, min, n}] ; 

Clear [rhs] ; 

« "aOl.txt"; 

<< "vOl.txt"; 
k= V 1+ 2a[0] ; 

2i+l 

Table[u[i, t] = ^ B[i, j] Cos[j t] , {i, 0, n, 2)] ; 

3*1 
A j=2 

2i+l 

Table[u[i, t] = ^ B[i, j] Sin[j t] , {i, 1, n, 2)] ; 

3=1 
A js2 

B« "bfile.txt"; 
w<< "wfile.txt"; 
b<< "crit.txt"; 

The next section calculates the four new terms, and are just like the equivalent sections of 
starter.nb except for using MyExpand[ ] instead of Expand[ ]: 

i= min; 

Print [ "VJbrking on term: ”, i] 

tenp = i:pdate[tenp, i] ; 

rhs = M^!:q>and[tenp|[i- min+ 1]] ; 
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si = Coefficient [rhs, Sin[k t] ] ; 

cl = Coefficient [rhs, Cos[k t] ] ; 

sol = Solve[ {si == 0, cl == 0 } , (a[i] , s[i - 1] ) ] ; 

a[i] = a[i] /. sol|[l]i ; 

s[i- 1] = s[i- 1] /. sol[l]| ; 

rhs = rhs / . sol|[l]i ; 

v[i, t] = rnYsolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] / . t-^ 0) == 1, c[i] ] [Ij ; 
i = min+ 1; 

Print ["Vforking on term: ", i] 

tenp= \jpdat/a[teap, i] ; 

rhs = I4^iq>and[tenp|[i- min+ Ij] ; 

si = Coefficient [rhs, Sin[k t] ] ; 

cl = Coefficient [rhs, Cos[k t] ] ; 

sol= Solve[(sl == 0, cl == 0} , {a[i] , s[i - 1]}] ; 

a[i] = a[i] /. sol[l]i ; 

s[i-l] = s[i-l] / . sollll ; 

rhs = rhs / . sollll ; 

v[i, t] = rrysolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[(getV[i] / . t-^ 0) == 1, c[i]][[l]l; 
i = min+ 2; 

Print [ "Vforking on term: ", i] 

tenp= V5x3ate[tenp, i] ; 

rhs = M^:q>and[ tenpli - min + IJ ] ; 

si = Coefficient [rhs, Sin[k t] ] ; 

cl = Coefficient [rhs, Cos[k t] ] ; 

sol= Solve[{sl== 0, cl== 0}, (a[i] , s[i- 1]}] ; 

a[i] = a[i] /. solllj ; 

s[i- 1] = s[i - 1] / . solllj ; 

rhs = rhs / . solllj ; 

v[i, t] = nysolve[i, rhs] ; 

v[i, t] = v[i, t] / . Solve[ <^tV[i] / . t-^ 0) == 1, c[i] ] [Ij ; 
i = min+ 3; 

Print [ "Wbrking on term: ", i] 

tenp= i53date[tenp, i] ; 

rhs = M^I:qpand[tenpli- min+ Ij] ; 

si = Coefficient [rhs, Sin[k t] ] ; 

cl = Coefficient [rhs, Cos[k t] ] ; 

sol = Solve[{sl == 0, cl == 0} , {a[i] , s[i - 1] } ] ; 
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a[i] = a[i] /. sol|[l]| ; 
s[i- 1] = s[i- 1] /. solilj ; 
rhs = rhs / . solflj ; 
v[i, t] = niYSolve[i, rhs] ; 

v[i, t] = v[i, t] /. Solve[ (getV[i] / . t-> 0) == 1, c[i]]|[l]i; 



Finally save the results as above: 



A= getA[n] ; 

V[t] = getV[n] ; 
DeleteFile [ "aOl . txt" ] 
Save[**a01-txt", a] 
DeleteFile [ "vOl . txt" ] 
Save [ "vOl . txt" , v] 
DeleteFile [ "nOl . txt" ] 
Save[ "nOl . txt" , n] 



APPENDIX C: SOLUTION TO VDP EQUATION 



The analytical solution to Equation (1-1) through orders’”: 

x[t] =2 Cos[t] + /_ Cos[t] ^ ^ Cos[3 1] — — COs[5 1] ] + 

^ 8 16 96 ’ 

4/73Cos[t] 47Cbs[3t] 1085Cbs[5t] 2149Cbs[7t] 6lCos[9t]\ 

^ ' 12288 1536 27648 110592 20480 ' 

6 / 6479 Cos [ t] 48437 Cos [ 3 1] 259945 Cbs [ 5 1] 

^ ' 663552^'" 35389440 31850496 

4253767 Cos [7 t] 480523 Cbs [9 1] 4937537 Cos [11 1] 715247 Cbs [13 1] \ 

398131200 73728000 2654208000 3715891200 ^ 

g / 377080601 Cos [ t] 36921629 Cos [ 3 1] 4094345 Cos [ 5 1] 

^ 1712282664960 71345111040 9172942848 

107409503789 Cos [ 7 1] 188691979247 Cos [9 1] 2067847855711 Cbs [11 1] 

45864714240000 59454259200000 936404582400000 

1114925731231 Cos [ 13 1] 784260289 Cos [ 15 1] 392636471 Cos [ 17 1] ^ 

1310966415360000 4661213921280 "" 29964946636800 ’ 

^10 / 44898976356847 Cos [t] 4475878408049 Cbs [3 1] 

^ ' 3020466620989440000 50341110349824000 

511229379671 Cos [ 5 1] 21294508407929 Cos [ 7 1] 1838562603094127 Cos [ 9 1] 

2958824445050880 165112971264000000 2621932830720000000 

1299787042060760957 Cbs [ 11 1] 11188471201628174701 Cos [ 13 1] 

1321454146682880000000 14800286442848256000000 

27736422934204291 Cos [ 15 1] 148525877527816577 Cos [ 17 1] 

78934861028524032000 1522315176978677760000 

1600597132693073 Cos [ 19 1] 29654422883 Cbs [21 1] ^ 

108736798355619840000 32288758824960000 ’ 

6 ( — — I5J .. - — Sin[3t]) +6^ Sin[3t] - — Sin[5t] + - 2 — Sin[7t]) + 

'4 4 '' 256 256 57 6 576 > 

5, 12971 Sin[t] 2591Sin[3t] 

^ ' 4423680 294912 

52885Sin[5t] 110621 Sin[ 7 1] 7457Sin[9t] 5533 Sin[ lit] ^ 

2654208 6635520 "" 1228800 7372800 > 

7 / 33114653 Sin[t] 82937Sin[3t] 1939625 Sin[ 5 1] 1061235889 Sin[ 7 1] 

^ ' 44590694400 212336640 764411904 191102976000 

179467921 Sin[ 9 1] 223342933 Sin[ 11 1] 195050323 Sin[ 13 1] 138697 Sin[ 15 1] ^ 

35389440000 "" 92897280000 346816512000 * 2774532096 ^ ^ 

9 r 9822191352131 Sin[t] 89866136849 Sin[ 3 1] 1702014109 Sin[ 5 1] 

^ ' 251705551749120000 342456532992000 12328435187712 

2088126932941 Sin[ 7 1] 41070940258801 Sin[ 9 1] 1694102483205143 Sin[ 11 1] 

2751882854400000 24970788864000000 1048773132288000000 

662514241768639 Sin[ 13 1] 507370714969 Sin[ 15 1] 15080405867887 Sin[ 17 1] 

734141192601600000 1740186530611200 "" 302046662098944000 

466445839 Sin[ 19 1] ^ 



134842259865600 
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Frequency Expansion through order 20: 

^ 1 , 35 e® 678899 

"l6 ^ 3072 * 884736 ' 5096079360 ^ 

28160413 e^O 16729607288111 5722795344767278507 

2293235712000 3698530556313600000 ” 5219366321069752320000000 ' 

1846779765852173498887007 3202506386889338212533493973243 

14731139504587268947968000000000 ^ 41577168137747107878744883200000000000 
4038372199485064617691118289043994731 e20 
3872464178615255430139995425341440000000000000 



NUMERICAL RESULTS 

These lists of numbers represent the coefficients of the cosine or sines of the VDP 
solution. c[2] for example shows that X 2 (t) = -0.125 Cos(t) + 0.188 Cos(t) -0.0521 
Cos(t). 

Cosine Coefficients through order 25: 

ctO] 

{ 2 .} 

c[2] 

{-0.125, 0.188, -0.0521} 
c[4] 

(0.00594, -0.0306, 0.0392, -0.0194, 0.00298} 

c[6] 

(0.000976, 0.00137, -0.00816, 0.0107, -0.00652, 0.00186, -0.000192} 

c[8] 

(-0.00022, 0.000518, 0.000446, -0.00234, 0.00317, -0.00221, 0.00085, -0.000168, 
0.0000131} 

c[10] 

(-0.0000149, -0.0000889, 0.000173, 0.000129, -0.000701, 0.000984, -0.000756, 

0.000351, -0.0000976, 0.0000147, -9.18x10"'^} 

c[12] 

(0.000011, -0.0000129, -0.000034, 0.000059, 0.0000364, -0.000214, 0.000313, 

-0.00026, 0.000139, -0.0000479, 0.0000104, -1.26x10'®, 6.56x10'®} 

c[14] 

(-3.84x10'’', 5.51x10'®, -4.1x10'®, -0.0000121, 0.0000203, 9.83x10'®, -0.0000661, 
0.000101, -0.00009, 0.0000533, -0.0000216, 5.91x10'®, -1.04x10'®, 1.06x10''', 
-4.75x10'®} 

c[16] 
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{-5.56x10 , 1.1x10 2.13x10"^, -1.4x10-^ -4.24x10'^ 7.02x 10 ^ 



2.46x 


10'®, 


-0.0000206, 


0.000033, -0 


.0000312, 


, 0 


.0000201, ■ 


-9.22x 


10'®, 3. 


xlO ®, 


-6.79 


xlO'® 


, 1.01x10®, 


-8.88x10 ®, 


3.47x10- 


-lOj 










c[18] 


10'®, 


-3.22x10'®, 
















(8.17x 


-8.11x10'®, 


7.85x10“ 


■ 7 

f 


-4.79x 10' 


\ -1.46x10'^, 




2.43x 


10'®, 


5.36x10'®, - 


■6.42x10'®, 0 


.0000108 




0.0000108, 


7.52x 


10'®, 




-3.81 


xlO'® 


, 1.42x10'®, 


-3.85x10'®, 


7.42x10" 


• 8 

f 


-9.55x 10'' 


^ 7.35 


X lO''^®, 


-2.55x10 


c[20] 






•1.22x10 ®, - 














(2.46x 


10'®, 


2.84x10'®, - 


1.12x10" 


• 8 

r 


2.85x10'®, 


-1.65 


X 10'®, - 


4.96x10'' 


8.45x 


10'®, 


8.07x10'®, - 


■2.01x10'®, 3 


.57x10"^ 


t ~ 


3.76x10'®, 


2.78x 


10'®, -1 


.53x10'®, 


6.37x 


10'®, 


-2. xlO ®, 4, 


,65x10'®, -7. 


78xl0-^ 


8. 


CO 

K) 

X 

o 


-6.05x 


10'®®, 1, 


.89x10®^ 



c[22] 

{-7.83x10'®, 1.68x10'®, 1.41x10'®, -4.49x10 ®, -6.79x10'®, 1.02x10'^ -5.82x10'®, 
-1.67x10'^ 2.94x10’’^, -7.04x10'®, -6.26xl0'\ 1.18xl0■^ -1.3x10'®, 

1.02x10'®, -6.05x10'^ 2.76xl0'\ -9.72x10'®, 2.63x10'®, -5.35x10'®, 

7.91x10'^®, -8.01x10'^^, 4.95x10'^, -1.41x10'^^} 

c(24] 

{-7.03x10'^®, -3.68x10'®, 6.03x10'®, 5.84x10'®, -1.63x10'®, -3.2x10'®, 

3.61x10'®, -2.08x10'®, -5.57x10'®, 1.02x10''', -1.41x10'®, -1.94x10'^ 

3.89x10'®, -4.51x10'®, 3.74x10 ®, -2.36x10'®, 1.16x10'®, -4.52x10'®, 

1.38x10'®, -3.28x10'®, 5.92x10'^®, -7.84x10'^®, 7.17x10'®^, -4.03x10'®®, 1.05x10'®^} 



Sine Coefficients through order 25: 

s(l] 

{0.75, -0.25} 
s[3] 

(-0.0273, 0.082, -0.0608, 0.0122} 

s[5] 

{-0.00293, -0.00879, 0.0199, -0.0167, 0.00607, -0.00075} 

s(7] 

{0.000743, -0.000391, -0.00254, 0.00555, -0.00507, 0.0024, -0.000562, 0.00005} 

s[9] 

{0.000039, 0.000262, -0.000138, -0.000759, 0.00164, -0.00162, 0.000902, 
-0.000292, 0.0000499, -3.46x10'®} 

S[ll] 

{-0.000035, -9.57x10'®, 0.0000953, -0.0000478, -0.000227, 0.000505, -0.000527, 
0.000332, -0.000132, 0.000032, -4.32x10'®, 2.45x10'®} 

s[13] 

{1.51x10'®, -0.0000106, -4.09x10'®, 0.0000335, -0.0000174, -0.0000683, 0.000158, 
-0.000175, 0.000121, -0.0000556, 0.000017, -3.3x 10'®, 3.67x10'®, -1.76x10'®} 

s[15] 
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{1.72x10'®, 1.56x10'®, -4.06x10'®, -1.73x10'®, 0.0000115, -6.5x10'®, 

-0.0000206, 0.0000501, -0.0000583, 0.0000437, -0.0000225, 8.14x10'®, 

-2.02x10 ®, 3.26x10'’', -3.08x10'®, 1.28x10'^} 

sfl7] 

{-2.66x10'’', 4.13x10'’', 6.46x10'’', -1.46x10'®, -6.44x10''', 3.94x10'®, 

-2.44x10'®, -6.16x10'®, 0.000016, -0.0000196, 0.0000157, -8.91x10'®, 

3.66x10'®, -1.09x10'®, 2.26x10'’', -3.12x10'®, 2.56x10'^ -9.4x10'^^} 

s[19] 

{-7.43x10'®, -1.34x10'’', 1.57x10'’', 2.59x10'^, -5.16x10'^ -2.24x10'’', 1.35x10'®, 
-9.17xl0'\ -1.83x10'®, 5.15x10'®, -6.62x10'®, 5.62x10'®, -3.45x10'®, 

1.58x10'®, -5.37x10'”', 1.35x10'”', -2.41x10 ®, 2.91x10'^ -2.11x10'^°, 6.94x10'^} 

s[21] 

{2.47x10'®, -1.02x10'®, -5.54x10'®, 5.45x10'®, 9.78x10'®, -1.81x10'”', 

-7.49x10'®, 4.59x10'”', -3.43x10'”', -5.37x10'”', 1.66x10'®, -2.24x10'®, 

2.01x10'®, -1.32x10'®, 6.59x10'”', -2.51x10'”', 7.3x10'®, -1.59x10'®, 

2.49x10'®, -2.66x10'^®, 1.73x10'’-^, -5.15x10'^} 

s[23] 

{2.02x10'®, 9.52x10'®, -3.45x10'®, -2.18x10'®, 1.86x10'®, 3.59x10'®, 

-6.37x10'®, -2.4x10'®, 1.56x10'”', -1.28x10'^, -1.54x10'^, 5.34x10'”', 

-7.57x10'”', 7.15x10'”', -5.X10'”', 2.69x10'”', -1.13x10'”', 3.69x10'®, 

-9.34x10'®, 1.79x10'®, -2.5x10'^®, 2.4x10'^, -1.41x10'''^, 3.84x10'^^} 

s[25] 

{-1.86x10'®, -4.12x10'^®, 3.91x10'®, -8.87x10'^®, -8.2x10'®, 6.36x10'®, 1.29x10'®, 
-2.23x10'®, -7.37x10'®, 5.29x10'®, -4.73x10'®, -4.31x10'®, 1.72x10'”', 

-2.56x10'”', 2.54x10'”', -1.87x10'”', 1.08x10'”', -4.9x10'®, 1.77x10'®, 

-5.08x10'®, 1.14x10'®, -1.94x10'^®, 2.45x10'^, -2.14x10'^, 1.15x10'^^, -2.88x10'^} 



APPENDIX D ; SOLUTIONS TO THE VARIATIONAL EQUATION 



This appendix includes the solutions to the variational equation through 8^ for the cases 
covered in Chapter 5. Note that these solutions can be used to check if the code in 
Appendix C is working correctly. 

Curves out of Ao=0: 
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Curves out of Ao=3/2: 
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Curves out of Ao=4: 
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359 3 . \ / 

Sn[3 1] SiiilS t] + - 

240 16 V V 



329Cos[t] 

640 



3Cos[t] 3 3 3Sm[t] 359 

+ — Cos[3 1 ] + — Cos[ 5 1 1 H + 

8 16 16 8 

233Cos[3t] 479Cos[5t] 27 51Sin[t] 73913Sin[3t] 15 27 \ 

• + 77:7 Cos[7 1] + — — — + ^“(3 1] - -777- Sin[7 1] I + 



1 - 



1280 
89813 Cos[t] 



1280 



640 



76151 Cos[3t] 139 

+ H Cos[5tJ + 

230400 600 



128 57600 

67Cos[7t] 



256 



640 



153600 230400 600 2048 

469Cos[9t] 831Sin[t] 82393Sin[3t] 679Sin[5t] 5311Sm[7t] 469Sin[9t]\ 

+ + + + 1 + 



46080 



2048 



57600 



5120 



51200 



46080 



1 - 



385217 Cos[t] 105987247 Cos[3 1] 

+ + 



614400 309657600 

389277 Cos[5t] 30401 Cos[7t] 327931 Cos[9t] 547Cos[llt] 273059Sin[t] 



1638400 


409600 


11059200 


215040 


614400 


37616912581 Sin[3t] 


18809 Sin[5t] 


249571 Sm[7t] 


27691 Sin[9t] 


547Sin[llt] ^ 


23224320000 


327680 


6144000 


2211840 


215040 ) 


42447078481 Cos[t] 


33657149027 Cos[3 1] 9222921259 Cos[5 1] 

• + + + 


1330129 Cos[7t] 



61931520000 92897280000 30%5760000 39321600 

46069577 Cos[9t] 88301 Cos[l It] 89371 Cos[13t] 400303339 Sin[t] 82024682603 Sin[ 3 1] 



9289728000 20643840 137625600 825753600 46448640000 

13130773 Sin[5t] 16057111 Sm[7t] 1 1346613 Sin[ 9 1] 883481 Sin[ 11 1] 89371 Sin[13t] ^ 

206438400 589824000 371589120 103219200 137625600 ) 



APPENDIX E : STABILITY TRANSITION CURVES 



This appendix includes the stability transition curves found in Chapter 5. Curves are given 
analytically through order 10 and numerically through order 30. 

Curves out of Ao=0; 

e le* 3e5 247 11657€^ 122481 

LoweieA = - — + — + — + + 

2 8 32 384 1024 442368 21233664 2548039680 

65987 129 248740367 

611529523200 36691771392000 

Lower sA = 

-0.5e + 0. 125e^ + 0.03125€^ - 0.018229€'' + 0.0029297€^ - 
0,00055836e^ - 0,00054899e'^ + 0.00048069 - 0.00010791 + 6.7792 x 10"^ e’® + 

0 . 000023325 - 0.000020161 + 3 . 8825 X 10 "^ e'^ + 7 . 1301 x lO'^e*'*- 1 . 1828 x 10 "^ + 

8 . 9935 X 10 -’e‘^- 1 . 1027 x lO-’e'"^ - 9 . 9766 x 10 '* + 6 . 0926 x lO'^e*^ - 3 . 6467 x 10 "*e 2 '^ + 
1 . 5256 X 10 '‘U^‘ + 8 . 7439 X 10 '^ e“ - 2 . 8529 x 10 '^ + 9 . 8924 x 10 '‘°€^ + 3.551 x 10 ' - 
6.2361 X 10 -‘°e 26 + 9 . 9193 x 10 " + 2 . 5033 x 10 "“ 3 . 8338 x 10 '“ + 3 . 6974 x 10 '" 
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Upper sA = 

0.5 £ + 0. 125 - 0.03 125 - 0.018229e‘» - 0.0029297 - 
0.00055836eO + 0.000548996'^ + 0.000480696* + 0.00010791 6^ + 6.7792x lO'^ £l0 - 
0.0000233256” - 0.000020161 €>2 -3.8825 x 10-^6” + 7.1301 x 10-'^£l'i + j ig2gx 10*^ 6*5 + 
8.9935X 10-2616+ 1.1027X 10-2612 - 9. 9766 x 1Q-* c** - 6.0926 x 10-*6l^ - 3.6467 x 10-*e20- 
1.5256X 10-”e2i 8 7439X 10-^622 + 2.8529x 10-^ £23 + 9.8924x 10-lle24 - 3.551 x 10-1H625 - 

6.2361 X 10-1“ £26- 9.9193 X iQ-W^n +2.5033x IQ-” £2* + 3.8338 x 10-”629 + 3.6974x lO’” 6^“ 



Curves out of Ao=3/2: 
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Curves out of Ao=4: 

5e2 2877 s'* 2537e5 54707 11663959e2 

40960 737280 ~ 15728640 14155776000 ~ 

43486217423 e9 28041629%941e‘" 

1223059046400000 14611478740992000 

Lower sA = 
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0.00082397£2-0.0013609e^ -0.000035555e9 + 0.00019192e*”- 0.000016118e“ + 
0.000037433e*2+ 1,4085 X 10-^ e“ - 0,00001 1912e*'‘ + 2.053 x IO’^e*^ - 6.2434 x IO-'^e*^- 
6.8113X 10-^e“+ 7,0367 X IO-’e** + 1.5598x IO-^e*^- 4.6186x IO-^e^” + 2.0464 x IO-^e^* - 
3.6947X 10'^e22-2.1736x 10’^ e^^ +7.4707x 10 "^e 2^+ 9.323x 10 ""e^^ + 1.5107x 10 ‘^e26 + 
1.8088X 10-’”e22 - 6.9398X IO-^^e^* - 2.4687x 10-“ ^ - 2.0128x lQ-“ 
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4. - 0.28125e2 + 0,078125e^ + 0.070239e'‘ - 0.003441 e^ - 0.0034782e6 - 
0.00082397 e’^ - 0.0013609e^ + 0,000035555 e^ + 0.00019192e*® + 0.000016118e“ + 

0.000037433 e“ - 1.4085 x 10-6 e“ _ o.000011912e 1'* - 2.053 x 10 -'^e“ - 6.2434x 10-'^el6 + 

6.81 13x 10 -^e“+ 7.0367 X IO-'^e*^ - 1.5598x 10 -^e> 9 - 4.6186x IO-^e^O- 2,0464 x 10-9 e^i _ 

3.6947 X 10-^e“+ 2.1736X 10-9 e^ + 7.4707 x \ q -9 ^ 10- “e^^ + 1.5 107 x 10-9 e26_ 

1.8088X 10-16 e2'7 - 6.9398X IO-^^e^* + 2.4687x 10-“ e^ - 2.0128x 10-“ e^” 
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